/* Crank_Nicolson.h
 *
 *	Copyright (C) 2008 CRIMERE
 *	Copyright (C) 2008 Jean-Marc Mercier
 *	
 *	This file is part of OTS (optimally transported schemes), an open-source library
 *	dedicated to scientific computing. http://code.google.com/p/optimally-transported-schemes/
 *
 *	CRIMERE makes no representations about the suitability of this
 *	software for any purpose. It is provided "as is" without express or
 *	implied warranty.
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 *
 */
#include <src/math/norms.h>
#include <src\algorithm\Evolution_Algorithm\TimeStepAlgorithm.h>
#include <src\math\Operators\Delta.h>

#if !defined(_CrankNicolson_)
#define _CrankNicolson_


template <class T>
struct CrankNicolson_Interface
{
	typedef CrankNicolson_Interface								type;
	CrankNicolson_Interface() {};

	virtual OS_double CFL() = 0;
	virtual bool criteria() = 0;
	virtual void gstar() = 0;
	virtual void vstar() = 0;

};


/*! \brief Crank Nicolson scheme

Crank_Nicolson generic schemes solved by a fixed point algorithm.
Given the initial condition \f$S^{t}\f$, a time increment \f$\delta_t\f$, and a non linear operator $Op(S)$,
we define a Crank_Nicolson scheme as the solution of the following non linear equation 
\f[
S^{n+1} = S^{t} + \delta_t Op(t^{1/2} ,S^{n+1/2) )
\f]
where \f$S^{n+1/2) = \frac {S^{t)+S^{n+1}}}{2}\f$ and \f$t^{1/2} = t + \Delta_t /2\f$.
the fixed point algorithm try to solve the previous equation by the following algorithm

\f[ S^{n+1} = S^{t) + \delta_t Op(t^{1/2} ,S^{n+1/2) )\f]

\f[ S^{0,n+1} : = S^{t); \f]
for \f$ k>=1 \f$, loop
\f[S^{k,n+1/2) := \frac {S^{t)+S^{k-1,n}{2}\f] 
\f[S^{k,n+1} : = S^{t) + \Delta_t Op(t^{1/2} ,S^{k,n+1/2) )\f]

until a convergence holds (Predictor), default being the \ell^2 norm

\f[|S^{k+1,n+1}-S^{k,n+1}|_{\ell^2}^2 < \epsilon, \text{ where } |S|_{\ell^2}^2 = \frac{1}{2N} \sum_{i=1,..N} |S|_i ^2\f]

This algorithm always converges for \f$\Delta\f$ small enough, under reasonnable conditions over the Operator Op.
However, much care has to be taken for the choice of the CFL condition, and we embed a simple time adaptive technique :
We verify \Delta \le cfl at the beginning of the algorithm.
If the algorithm seems to diverge (more than 10 iterations), then we try to divide the step time \Delta by 2 and reset the algorithm.
If the algorithm resumes to the Forward Euler Algorithm (only one iteration), then the step time is multiplied by two.
*/
template <	class TimeTypeOperator >
class Crank_Nicolson : public TimeStepAlgorithm<TimeTypeOperator>
	,public CrankNicolson_Interface<TimeTypeOperator>
{
public:
	typedef typename TimeStepAlgorithm<time_operator_type>			Base;
	typedef typename Crank_Nicolson<time_operator_type>				algorithm_type;
	typedef typename TimeTypeOperator::argument_type				argument_type;
	typedef typename TimeTypeOperator::argument_base_type				argument_base_type;
	virtual ~Crank_Nicolson(){};
	Crank_Nicolson() : cfl_adjust_(1.), n_(0) {};
	virtual void gstar() {
		Op_iterate = Op_(t_12,S_k12);
	};
	virtual void vstar() {
		S_k12   = (*get_state_variable()+iterate)/2.;
	};

	result_type operator()(const time_type & from, time_type & to, argument_type S_k) {

		current_time_ = from;
		nb_iterations_ = 0.;
		delta_t_	= to-from ;
		iterate		= S_k;
		S_k12		= S_k;
		t_12        = from;
		Op_iterate  = Op_(t_12,S_k12);
		delta_t_ =  std::min(delta_t_, cfl_adjust_*CFL());
//		OS_DYNAMIC_CHECK(delta_t_ > THRESHOLD,"Crank_Nicolson : the time step is too low. Check the CFL function.");
		vstar();
		gstar();
		iterate    = S_k+delta_t_*Op_iterate;
		t_12		= from + delta_t_ /2. ;
		bool Convergent = false;
		do {
start:				
				nb_iterations_++;

				Op_S_k1 = Op_iterate;
				S_k1    = iterate;
				vstar();
				gstar();
				iterate = S_k + delta_t_*Op_iterate;

				if ( criteria() ) {
					break;
				}
				if (nb_iterations_ > 10 ) {
					adjust_cfl(from,S_k);
				}
		} while (true);
		if (nb_iterations_ < 3) {
//From one side, we should never overpass the cfl condition, that ensures stability.
			if (cfl_adjust_ < 1)		
				cfl_adjust_ *= 2.; // multiply CFL constant
//From other side, it may be really hard to determine a correct CFL function. 
// May be a external parameter defaulted to false is a good choice.

		}
end:
		to = from + delta_t_;
		++n_;
		return iterate;
	}
	static std::string get_name(){return " Crank Nicolson delta_t u = ";};
protected :
	virtual void adjust_cfl(const time_type & from, argument_type S_k){
		OS_DYNAMIC_CHECK(cfl_adjust_ > THRESHOLD,"Crank_Nicolson : problem in Crank Nicolson scheme");
		cfl_adjust_ /= 2.; // divide CFL constant
		delta_t_ /= 2.; 
		t_12		= from + delta_t_ / 2.;
		iterate = *get_state_variable();
		vstar();
		gstar();
		iterate =  (S_k+delta_t_* Op_iterate);
		nb_iterations_ = 0;
	}
	indice_type n_;
	OS_double cfl_adjust_;
	time_type t_12, delta_t_;
	time_type current_time_;
	argument_base_type S_k12, iterate, Op_iterate, S_k1, Op_S_k1;
	std::size_t nb_iterations_;
};


#endif
